High-Order Discontinuous Galerkin Finite Element Methods with Globally Divergence-Free Constrained Transport for Ideal MHD

The modification of the celebrated Yee scheme from Maxwell equations to magnetohydrodynamics is often referred to as the constrained transport approach. Constrained transport can be viewed as a sort of predictor-corrector method for updating the magnetic field, where a magnetic field value is first predicted by a method that does not preserve the divergence-free condition on the magnetic field, followed by a correction step that aims to control these divergence errors. This strategy has been successfully used in conjunction with a variety of shock-capturing methods including WENO, central, and wave propagation schemes. In this work we show how to extend the basic CT framework to the discontinuous Galerkin finite element method on both 2D and 3D Cartesian grids. We first review the entropy-stability theory for semi-discrete DG discretizations of ideal MHD, which rigorously establishes the need for a magnetic field that satisfies the following conditions: (1) the divergence of the magnetic field is zero on each element, and (2) the normal components of the magnetic field are continuous across element edges/faces. In order to achieve such a globally divergence-free magnetic field, we introduce a novel CT scheme that is based on two ingredients: (1) we introduce an element-centered magnetic vector potential that is updated via a discontinuous Galerkin scheme on the induction equation; and (2) we define a mapping that takes element-centered magnetic field values and element-centered magnetic vector potential values and creates on each edge /face a representation of the normal component of the magnetic field; this representation is then mapped back to the elements to create a globally divergence-free element-centered representation of the magnetic field. For problems with shock waves, we make use of so-called moment-based limiters to control oscillations in the conserved quantities.Comment: 26 pages, 6 figure

A Class of Quadrature-Based Moment-Closure Methods with Application to the Vlasov-Poisson-Fokker-Planck System in the High-Field Limit

Quadrature-based moment-closure methods are a class of approximations that replace high-dimensional kinetic descriptions with lower-dimensional fluid models. In this work we investigate some of the properties of a sub-class of these methods based on bi-delta, bi-Gaussian, and bi-B-spline representations. We develop a high-order discontinuous Galerkin (DG) scheme to solve the resulting fluid systems. Finally, via this high-order DG scheme and Strang operator splitting to handle the collision term, we simulate the fluid-closure models in the context of the Vlasov-Poisson-Fokker-Planck system in the high-field limit. We demonstrate numerically that the proposed scheme is asymptotic-preserving in the high-field limit.Comment: 24 pages, 5 figure

Outflow Positivity Limiting for Hyperbolic Conservation Laws. Part I: Framework and Recipe

Numerical methods for hyperbolic conservation laws are needed that efficiently mimic the constraints satisfied by exact solutions, including material conservation and positivity, while also maintaining high-order accuracy and numerical stability. Discontinuous Galerkin (DG) and WENO schemes allow efficient high-order accuracy while maintaining conservation. Positivity limiters developed by Zhang and Shu ensure a minimum time step for which positivity of cell average quantities is maintained without sacrificing conservation or formal accuracy; this is achieved by linearly damping the deviation from the cell average just enough to enforce a cell positivity condition that requires positivity at boundary nodes and strategically chosen interior points. We assume that the set of positive states is convex; it follows that positivity is equivalent to scalar positivity of a collection of affine functionals. Based on this observation, we generalize the method of Zhang and Shu to a framework that we call outflow positivity limiting: First, enforce positivity at boundary nodes. If wave speed desingularization is needed, cap wave speeds at physically justified maxima by using remapped states to calculate fluxes. Second, apply linear damping again to cap the boundary average of all positivity functionals at the maximum possible (relative to the cell average) for a scalar-valued representation positive in each mesh cell. This be done by enforcing positivity of the retentional, an affine combination of the cell average and the boundary average, in the same way that Zhang and Shu would enforce positivity at a single point (and with similar computational expense). Third, limit the time step so that cell outflow is less than the initial cell content. This framework guarantees essentially the same positivity-preserving time step as is guaranteed if positivity is enforced at every point in the mesh cell.Comment: 32 pages, 15 figure

The Regionally-Implicit Discontinuous Galerkin Method: Improving the Stability of DG-FEM

Discontinuous Galerkin (DG) methods for hyperbolic partial differential equations (PDEs) with explicit time-stepping schemes, such as strong stability-preserving Runge-Kutta (SSP-RK), suffer from time-step restrictions that are significantly worse than what a simple Courant-Friedrichs-Lewy (CFL) argument requires. In particular, the maximum stable time-step scales inversely with the highest degree in the DG polynomial approximation space and becomes progressively smaller with each added spatial dimension. In this work we introduce a novel approach that we have dubbed the regionally implicit discontinuous Galerkin (RIDG) method to overcome these small time-step restrictions. The RIDG method is based on an extension of the Lax-Wendroff DG (LxW-DG) method, which previously had been shown to be equivalent to a predictor-corrector approach, where the predictor is a locally implicit spacetime method (i.e., the predictor is something like a block-Jacobi update for a fully implicit spacetime DG method). The corrector is an explicit method that uses the spacetime reconstructed solution from the predictor step. In this work we modify the predictor to include not just local information, but also neighboring information. With this modification we show that the stability is greatly enhanced; in particular, we show that we are able to remove the polynomial degree dependence of the maximum time-step and show how this extends to multiple spatial dimensions. A semi-analytic von Neumann analysis is presented to theoretically justify the stability claims. Convergence and efficiency studies for linear and nonlinear problems in multiple dimensions are accomplished using a MATLAB code that can be freely downloaded.Comment: 26 pages, 4 figures, 8 table

Ten-moment two-fluid plasma model agrees well with PIC/Vlasov in GEM problem

We simulate magnetic reconnection in the GEM problem using a two-fluid model with 10 moments for the electron fluid as well as the proton fluid. We show that use of 10 moments for electrons gives good qualitative agreement with the the electron pressure tensor components in published kinetic simulations.Comment: 13th International Conference on Hyperbolic Problems, June 201

Simulation of Fast Magnetic Reconnection using a Two-Fluid Model of Collisionless Pair Plasma without Anomalous Resistivity

For the first time to our knowledge, we demonstrate fast magnetic reconnection near a magnetic null point in a fluid model of collisionless pair plasma without resorting to the contrivance of anomalous resistivity. In particular, we demonstrate that fast reconnection occurs in an anisotropic adiabatic two-fluid model of collisionless pair plasma with relaxation toward isotropy for a broad range of isotropization rates. For very rapid isotropization we see fast reconnection, but instabilities eventually arise that cause numerical error and cast doubt on the simulated behavior

A high-order unstaggered constrained transport method for the 3D ideal magnetohydrodynamic equations based on the method of lines

Numerical methods for solving the ideal magnetohydrodynamic (MHD) equations in more than one space dimension must confront the challenge of controlling errors in the discrete divergence of the magnetic field. One approach that has been shown successful in stabilizing MHD calculations are constrained transport (CT) schemes. CT schemes can be viewed as predictor-corrector methods for updating the magnetic field, where a magnetic field value is first predicted by a method that does not exactly preserve the divergence-free condition on the magnetic field, followed by a correction step that aims to control these divergence errors. In Helzel et al. (2011) the authors presented an unstaggered constrained transport method for the MHD equations on 3D Cartesian grids. In this work we generalize the method of Helzel et al. (2011) in three important ways: (1) we remove the need for operator splitting by switching to an appropriate method of lines discretization and coupling this with a non-conservative finite volume method for the magnetic vector potential equation, (2) we increase the spatial and temporal order of accuracy of the entire method to third order, and (3) we develop the method so that it is applicable on both Cartesian and logically rectangular mapped grids. The evolution equation for the magnetic vector potential is solved using a non-conservative finite volume method. The curl of the magnetic potential is computed via a third-order accurate discrete operator that is derived from appropriate application of the divergence theorem and subsequent numerical quadrature on element faces. Special artificial resistivity limiters are used to control unphysical oscillations in the magnetic potential and field components across shocks. Test computations are shown that confirm third order accuracy for smooth test problems and high-resolution for test problems with shock waves.Comment: 29 pages, 5 figures, 3 table

A Simple and Effective High-Order Shock-Capturing Limiter for Discontinuous Galerkin Methods

The discontinuous Galerkin (DG) finite element method when applied to hyperbolic conservation laws requires the use of shock-capturing limiters in order to suppress unphysical oscillations near large solution gradients. In this work we develop a novel shock-capturing limiter that combines key ideas from the limiter of Barth and Jespersen [AIAA-89-0366 (1989)] and the maximum principle preserving (MPP) framework of Zhang and Shu [Proc. R. Soc. A, 467 (2011), pp. 2752--2776]. The limiting strategy is based on traversing the mesh element-by-element in order to (1) find local upper and lower bounds on user-defined variables by sampling these variables on neighboring elements, and (2) to then enforce these local bounds by minimally damping the high-order corrections. The main advantages of this limiting strategy is that it is simple to implement, effective at shock capturing, and retains high-order accuracy of the solution in smooth regimes. The resulting numerical scheme is applied to several standard numerical tests in both one and two-dimensions and on both Cartesian and unstructured grids. These tests are used as benchmarks to verify and assess the accuracy and robustness of the method.Comment: 20 page

Positivity-preserving discontinuous Galerkin methods with Lax-Wendroff time discretizations

This work introduces a single-stage, single-step method for the compressible Euler equations that is provably positivity-preserving and can be applied on both Cartesian and unstructured meshes. This method is the first case of a single-stage, single-step method that is simultaneously high-order, positivity-preserving, and operates on unstructured meshes. Time-stepping is accomplished via the Lax-Wendroff approach, which is also sometimes called the Cauchy-Kovalevskaya procedure, where temporal derivatives in a Taylor series in time are exchanged for spatial derivatives. The Lax-Wendroff discontinuous Galerkin (LxW-DG) method developed in this work is formulated so that it looks like a forward Euler update but with a high-order time-extrapolated flux. In particular, the numerical flux used in this work is a linear combination of a low-order positivity-preserving contribution and a high-order component that can be damped to enforce positivity of the cell averages for the density and pressure for each time step. In addition to this flux limiter, a moment limiter is applied that forces positivity of the solution at finitely many quadrature points within each cell. The combination of the flux limiter and the moment limiter guarantees positivity of the cell averages from one time-step to the next. Finally, a simple shock capturing limiter that uses the same basic technology as the moment limiter is introduced in order to obtain non-oscillatory results. The resulting scheme can be extended to arbitrary order without increasing the size of the effective stencil. We present numerical results in one and two space dimensions that demonstrate the robustness of the proposed scheme.Comment: 28 pages, 9 figure

Finite Difference Weighted Essentially Non-Oscillatory Schemes with Constrained Transport for Ideal Magnetohydrodynamics

In this work we develop a class of high-order finite difference weighted essentially non-oscillatory (FD-WENO) schemes for solving the ideal magnetohydrodynamic (MHD) equations in 2D and 3D. The philosophy of this work is to use efficient high-order WENO spatial discretizations with high-order strong stability-preserving Runge-Kutta (SSP-RK) time-stepping schemes. Numerical results have shown that with such methods we are able to resolve solution structures that are only visible at much higher grid resolutions with lower-order schemes. The key challenge in applying such methods to ideal MHD is to control divergence errors in the magnetic field. We achieve this by augmenting the base scheme with a novel high-order constrained transport approach that updates the magnetic vector potential. The predicted magnetic field from the base scheme is replaced by a divergence-free magnetic field that is obtained from the curl of this magnetic potential. The non-conservative weakly hyperbolic system that the magnetic vector potential satisfies is solved using a version of FD-WENO developed for Hamilton-Jacobi equations. The resulting numerical method is endowed with several important properties: (1) all quantities, including all components of the magnetic field and magnetic potential, are treated as point values on the same mesh (i.e., there is no mesh staggering); (2) both the spatial and temporal orders of accuracy are fourth-order; (3) no spatial integration or multidimensional reconstructions are needed in any step; and (4) special limiters in the magnetic vector potential update are used to control unphysical oscillations in the magnetic field. Several 2D and 3D numerical examples are presented to verify the order of accuracy on smooth test problems and to show high-resolution on test problems that involve shocks.Comment: 39 pages, 9 figures, 4 table